# DiD for placebo outcomes (sports / comedy viewing)

library(tidyverse)
library(magrittr)
library(data.table)
library(lubridate)
library(hms)
library(ggthemes)
library(stargazer)
theme_set(theme_few())

### SET WORKING DIRECTORY HERE ###
path_to_archive <- "replication/"
data_dir <- paste0(path_to_archive, "data/")
setwd(data_dir)
plot_dir <- paste0(path_to_archive, "plots/")
tables_dir <- paste0(path_to_archive, "tables/")

diff2 <- function(x) x[seq(1,length(x)-1,2)] - x[seq(2,length(x),2)]
inv_var_weight <- function(i,d) {
	weighted.mean(d[i, est], w=sqrt(d[i, n_t] * d[i, n_c]), na.rm=T)
}

total_size_weight <- function(i, d) {
	weighted.mean(d[i,est], w=d[i,n_t] + d[i,n_c], na.rm=T)
}

simple_avg <- function(i,d){
	mean(d[i,est], na.rm=T)
}

dd_wide <- readRDS(paste0(data_dir, "all_dd_data_sports.rds"))
s_cols <- grep("s_p", colnames(dd_wide), value=T)

# construct combined C group
all_controls <- dd_wide[group!="T", map(.SD, sum), .SDcols = s_cols, by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12)] %>%
	.[,group:="C"]

dd_wide %<>% 
	rbind(all_controls) %>%
	.[,n_c := n_c1 + n_c2 + n_c12] %>%
	.[group=="C", (s_cols) := map(.SD, ~ . / n_c), .SDcols = s_cols] %>% 
	.[group=="C1", (s_cols) := map(.SD, ~ . / n_c1), .SDcols = s_cols] %>% 
	.[group=="C2", (s_cols) := map(.SD, ~ . / n_c2), .SDcols = s_cols] %>% 
	.[group=="C1+C2", (s_cols) := map(.SD, ~ . / n_c12), .SDcols = s_cols] %>% 
	.[group=="T", (s_cols) := map(.SD, ~ . / n_t), .SDcols = s_cols] 


dd_long <- dd_wide %>% 
	gather(key="hour", val="view_mean", dplyr::starts_with("s_pre"), dplyr::starts_with("s_post")) %>%
	as.data.table

dd_long[, `:=`(pre_post = sub("s_(pre|post)_(\\d+)", "\\1", hour),
			   hour = sub("s_(pre|post)_(\\d+)", "\\2", hour) %>% as.numeric)]

dd_base <- 
	dd_long[group %in% c("C","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

# full sample point estimate, bootstrap
pt_est_full <- inv_var_weight(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_base), replace=T), dd_base))
results_full <- data.table(replicate = 0:1000, estimate = c(pt_est_full, bs_full))

pt_est_full_simple <- simple_avg(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full_simple <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_base), replace=T), dd_base))
results_full_simple <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple, bs_full_simple))

pt_est_full_size <- total_size_weight(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full_size <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_base), replace=T), dd_base))
results_full_size <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size, bs_full_size))

## separately for C1 / C2 / C3
dd_c1 <- 
	dd_long[group %in% c("C1","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c1 <- inv_var_weight(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_c1 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c1, bs_full_c1))

pt_est_full_simple_c1 <- simple_avg(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_simple_c1 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_simple_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c1, bs_full_simple_c1))

pt_est_full_size_c1 <- total_size_weight(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_size_c1 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_size_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c1, bs_full_size_c1))


dd_c2 <- 
	dd_long[group %in% c("C2","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c2 <- inv_var_weight(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_c2 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c2, bs_full_c2))

pt_est_full_simple_c2 <- simple_avg(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_simple_c2 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_simple_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c2, bs_full_simple_c2))

pt_est_full_size_c2 <- total_size_weight(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_size_c2 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_size_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c2, bs_full_size_c2))

dd_c12 <- 
	dd_long[group %in% c("C1+C2","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c12 <- inv_var_weight(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_c12 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c12, bs_full_c12))

pt_est_full_simple_c12 <- simple_avg(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_simple_c12 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_simple_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c12, bs_full_simple_c12))

pt_est_full_size_c12 <- total_size_weight(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_size_c12 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_size_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c12, bs_full_size_c12))

quantiles <- c(0.01,0.025,0.975,0.99)
ci_full_c0 <- results_full[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c1 <- results_full_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c2 <- results_full_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c12 <- results_full_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c0 <- results_full_size[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c1 <- results_full_size_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c2 <- results_full_size_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c12 <- results_full_size_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c0 <- results_full_simple[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c1 <- results_full_simple_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c2 <- results_full_simple_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c12 <- results_full_simple_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]

ci_full_split <- rbind(
	ci_full_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])
ci_full_split_size <- rbind(
	ci_full_size_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_size_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_size_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_size_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])
ci_full_split_simple <- rbind(
	ci_full_simple_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_simple_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_simple_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_simple_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])

### TABLE C.2
y <- rnorm(100)
x <- matrix(rnorm(100*nrow(ci_full_split)), nrow=100)
fake_model <- lm(y ~ x - 1)

dd_bootstrap_cis_split <- stargazer(list(fake_model, fake_model, fake_model),
		  title="Average Effects and Bootstrapped Confidence Intervals, Sports viewing placebo.",
		  style="apsr",
		  out=paste0(tables_dir, "dd_bootstrap_cis_sports.tex"),
		  label="tab:dd_boostrap_cis_sports",
		  coef=list(ci_full_split[,point_est], ci_full_split_size[,point_est], ci_full_split_simple[,point_est]), 
		  ci=T,
		  ci.custom=list(ci_full_split[, .(q_025, q_975)] %>% as.matrix,
		  				 ci_full_split_size[, .(q_025, q_975)] %>% as.matrix,
		  				 ci_full_split_simple[, .(q_025, q_975)] %>% as.matrix),
		  star.cutoffs = NA,
		  omit.stat = "all",
		  column.labels = c("Inv. Variance", "Total Devices", "Equally Weighted"),
		  dep.var.labels = "Effect (Sports Viewing Minutes)",
		  covariate.labels = ci_full_split[,split],
		  notes="\\parbox[t]{0.9\\linewidth}{Table reports average effects on sports program viewing (in minutes), weighted by 1) the geometric mean of treatment and control group size, 2) total devices in treatment and control groups and 3) a simple equally weighted average. Confidence intervals are the central 95\\% interval of 1000 bootstrap replicates within each subgroup. Group C1 is devices in control group which tuned out prior to ad time; C2 is devices in control group which tuned in after ad time; and C3 is devices in control group who tuned out before ad time, and then tuned back in after ad time. 'All' row includes all three types as controls.}",
		  notes.append=F
		  )

### Now for comedy ###


dd_wide <- readRDS(paste0(data_dir, "all_dd_data_comedy.rds"))
s_cols <- grep("s_p", colnames(dd_wide), value=T)

# construct combined C group
all_controls <- dd_wide[group!="T", map(.SD, sum), .SDcols = s_cols, by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12)] %>%
	.[,group:="C"]

dd_wide %<>% 
	rbind(all_controls) %>%
	.[,n_c := n_c1 + n_c2 + n_c12] %>%
	.[group=="C", (s_cols) := map(.SD, ~ . / n_c), .SDcols = s_cols] %>% 
	.[group=="C1", (s_cols) := map(.SD, ~ . / n_c1), .SDcols = s_cols] %>% 
	.[group=="C2", (s_cols) := map(.SD, ~ . / n_c2), .SDcols = s_cols] %>% 
	.[group=="C1+C2", (s_cols) := map(.SD, ~ . / n_c12), .SDcols = s_cols] %>% 
	.[group=="T", (s_cols) := map(.SD, ~ . / n_t), .SDcols = s_cols] 


dd_long <- dd_wide %>% 
	gather(key="hour", val="view_mean", dplyr::starts_with("s_pre"), dplyr::starts_with("s_post")) %>%
	as.data.table

dd_long[, `:=`(pre_post = sub("s_(pre|post)_(\\d+)", "\\1", hour),
			   hour = sub("s_(pre|post)_(\\d+)", "\\2", hour) %>% as.numeric)]

dd_base <- 
	dd_long[group %in% c("C","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

# full sample point estimate, bootstrap
pt_est_full <- inv_var_weight(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_base), replace=T), dd_base))
results_full <- data.table(replicate = 0:1000, estimate = c(pt_est_full, bs_full))

pt_est_full_simple <- simple_avg(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full_simple <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_base), replace=T), dd_base))
results_full_simple <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple, bs_full_simple))

pt_est_full_size <- total_size_weight(1:nrow(dd_base),dd_base)
set.seed(396843)
bs_full_size <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_base), replace=T), dd_base))
results_full_size <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size, bs_full_size))

## separately for C1 / C2 / C3
dd_c1 <- 
	dd_long[group %in% c("C1","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c1 <- inv_var_weight(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_c1 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c1, bs_full_c1))

pt_est_full_simple_c1 <- simple_avg(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_simple_c1 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_simple_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c1, bs_full_simple_c1))

pt_est_full_size_c1 <- total_size_weight(1:nrow(dd_c1),dd_c1)
set.seed(1850159)
bs_full_size_c1 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c1), replace=T), dd_c1))
results_full_size_c1 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c1, bs_full_size_c1))


dd_c2 <- 
	dd_long[group %in% c("C2","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c2 <- inv_var_weight(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_c2 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c2, bs_full_c2))

pt_est_full_simple_c2 <- simple_avg(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_simple_c2 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_simple_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c2, bs_full_simple_c2))

pt_est_full_size_c2 <- total_size_weight(1:nrow(dd_c2),dd_c2)
set.seed(1850159)
bs_full_size_c2 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c2), replace=T), dd_c2))
results_full_size_c2 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c2, bs_full_size_c2))

dd_c12 <- 
	dd_long[group %in% c("C1+C2","T"), .(news_min = sum(view_mean)/60), by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post)] %>%
	setkey(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c, group, pre_post) %>%
	.[, .(est = -diff2(diff2(.SD$news_min))) ,by = .(dma_code, channel, affiliate, program, ad_time_utc, n_t, n_c1, n_c2, n_c12,n_c)]

pt_est_full_c12 <- inv_var_weight(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_c12 <- map_dbl(1:1000, ~ inv_var_weight(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_c12, bs_full_c12))

pt_est_full_simple_c12 <- simple_avg(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_simple_c12 <- map_dbl(1:1000, ~ simple_avg(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_simple_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_simple_c12, bs_full_simple_c12))

pt_est_full_size_c12 <- total_size_weight(1:nrow(dd_c12),dd_c12)
set.seed(1850159)
bs_full_size_c12 <- map_dbl(1:1000, ~ total_size_weight(sample.int(nrow(dd_c12), replace=T), dd_c12))
results_full_size_c12 <- data.table(replicate = 0:1000, estimate = c(pt_est_full_size_c12, bs_full_size_c12))

quantiles <- c(0.01,0.025,0.975,0.99)
ci_full_c0 <- results_full[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c1 <- results_full_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c2 <- results_full_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_c12 <- results_full_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c0 <- results_full_size[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c1 <- results_full_size_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c2 <- results_full_size_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_size_c12 <- results_full_size_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c0 <- results_full_simple[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c1 <- results_full_simple_c1[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c2 <- results_full_simple_c2[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]
ci_full_simple_c12 <- results_full_simple_c12[,c(point_est=estimate[replicate==0], map(quantiles, ~ quantile(estimate[replicate!=0], probs=.)) %>% set_names(paste("q",c("01","025","975","99"), sep="_")))]

ci_full_split <- rbind(
	ci_full_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])
ci_full_split_size <- rbind(
	ci_full_size_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_size_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_size_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_size_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])
ci_full_split_simple <- rbind(
	ci_full_simple_c0[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="All")],
	ci_full_simple_c1[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C1")],
	ci_full_simple_c2[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C2")],
	ci_full_simple_c12[,.(point_est,q_01,q_025,q_975,q_99)][,`:=` (split="C3")])


### TABLE C.3
y <- rnorm(100)
x <- matrix(rnorm(100*nrow(ci_full_split)), nrow=100)
fake_model <- lm(y ~ x - 1)

dd_bootstrap_cis_split <- stargazer(list(fake_model, fake_model, fake_model),
		  title="Average Effects and Bootstrapped Confidence Intervals, Comedy viewing placebo.",
		  style="apsr",
		  out=paste0(tables_dir, "dd_bootstrap_cis_comedy.tex"),
		  label="tab:dd_boostrap_cis_comedy",
		  coef=list(ci_full_split[,point_est], ci_full_split_size[,point_est], ci_full_split_simple[,point_est]), 
		  ci=T,
		  ci.custom=list(ci_full_split[, .(q_025, q_975)] %>% as.matrix,
		  				 ci_full_split_size[, .(q_025, q_975)] %>% as.matrix,
		  				 ci_full_split_simple[, .(q_025, q_975)] %>% as.matrix),
		  star.cutoffs = NA,
		  omit.stat = "all",
		  column.labels = c("Inv. Variance", "Total Devices", "Equally Weighted"),
		  dep.var.labels = "Effect (Comedy Viewing Minutes)",
		  covariate.labels = ci_full_split[,split],
		  notes="\\parbox[t]{0.9\\linewidth}{Table reports average effects on comedy program viewing (in minutes), weighted by 1) the geometric mean of treatment and control group size, 2) total devices in treatment and control groups and 3) a simple equally weighted average. Confidence intervals are the central 95\\% interval of 1000 bootstrap replicates within each subgroup. Group C1 is devices in control group which tuned out prior to ad time; C2 is devices in control group which tuned in after ad time; and C3 is devices in control group who tuned out before ad time, and then tuned back in after ad time. 'All' row includes all three types as controls.}",
		  notes.append=F
		  )
